function [x fval] = emd(w1, w2, C)
%
% EMD   Earth Mover's Distance between two signatures
%    [X, FVAL] = EMD(w1, w2, C) is the Earth Mover's Distance
%    between two signatures with in W1 and W2. C is the cost matrix 
%    between two corresponding feature vectors.
%
%    EMD is formalized as linear programming problem in which the flow that
%    minimizes an overall cost function subject to a set of constraints is
%    computed. This implementation is based on "The Earth Mover's Distance
%    as a Metric for Image Retrieval", Y. Rubner, C. Tomasi and L. Guibas,
%    International Journal of Computer Vision, 40(2), pp. 99-121, 2000.
%
%    The outcome of EMD is the flow (X) which minimizes the cost function
%    and the value (FVAL) of this flow.
%
%    This file and its content belong to Ulas Yilmaz.
%    You are welcome to use it for non-commercial purposes, such as
%    student projects, research and personal interest. However,
%    you are not allowed to use it for commercial purposes, without
%    an explicit written and signed license agreement with Ulas Yilmaz.
%    Berlin University of Technology, Germany 2006.
%    http://www.cv.tu-berlin.de/~ulas/RaRF
%
%    modified by Shaobo Hou, 21/02/2010
%


% number of feature vectors
M = numel(w1);
N = numel(w2);

% inequality constraints
A1 = zeros(M, M, N);
A2 = zeros(N, M, N);
for m = 1:M
    A1(m, m, :) = 1;
end
for n = 1:N
    A2(n, :, n) = 1;
end
A1 = reshape(A1, [M M*N]);
A2 = reshape(A2, [N M*N]);
A = [A1; A2];
b = [w1(:); w2(:)];

% equality constraints
Aeq = ones(M+N, M*N);
beq = ones(M+N, 1) * min(sum(w1), sum(w2));

% lower bound
lb = zeros(1, M*N);

% linear programming
[x fval] = linprog(C, A, b, Aeq, beq, lb);
fval = fval / sum(x);
